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Abstract 

Valuing Guaranteed Lifelong Withdrawal Benefit (GLWB) has attracted significant attention from both 
the academic field and real world financial markets. As remarked by Forsyth and Vetzal [S] the Black 
and Scholes framework seems to be inappropriate for such long maturity products. They propose to use a 
regime switching model. Alternatively, we propose here to use a stochastic volatility model (Heston model) 
and a Black Scholes model with stochastic interest rate (Hull White model). For this purpose we present 
four numerical methods for pricing GLWB variables annuities: a hybrid tree-finite difference method and a 
hybrid Monte Garlo method, an ADI finite difference scheme, and a standard Monte Garlo method. These 
methods are used to determine the no-arbitrage fee for the most popular versions of the GLWB contract, 
and to calculate the Greeks used in hedging. Both constant withdrawal and optimal withdrawal (including 
lapsation) strategies are considered. Numerical results are presented which demonstrate the sensitivity of 
the no-arbitrage fee to economic, contractual and longevity assumptions. 
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1 Introduction 


In 2008 following the subprime crisis, financial markets have suffered the upheavals that have affected the 
entire world economy. Since then, these markets were extremely volatile: this situation could last a while, 
and perhaps become the new standard. After many failures, the gap between the different interest rate apply 
to different transmitters has become larger and larger, and a discussion on the identification of the risk-free 
rate is opened. The ECB and the Fed’s rate gradually declined, while rate on sovereign debt increased 
gradually. 

For customers, it is difficult to balance risk and return. In this context, clients seek protection for 
their savings, and the ability to take advantage of the positive changes in the market. With regard social 
problematic, following the increase in life expectancy, annuities for retirement dropped. 

The mission of insurance companies is to answer the request for protection and compensation of its 
customers. The solution is to provide the customer an investment account and cover its value with guarantees. 
These products are called Variable Annuities. In the words of Francois Robinet, CFO of AXA Life Invest, 
"These products, unit of account guaranteed will become a solution to solve the long-term investment 
problems with security, prepare for retirement". 

In this article, we consider a Guaranteed Lifelong Withdrawal Benefit (GLWB) annuity. We restrict 
our attention to a simplified form of a GLWB which is initiated by making a lump sum payment to an 
insurance company. This lump sum is then invested in risky assets, usually a mutual fund. The benefit 
base, or guarantee account balance, is initially set to the amount of the lump sum payment. The holder of 
the contract is entitled to withdraw a fixed fraction of the benefit base for life, even if the actual investment 
in the risky asset declines to zero. Upon the death of the contract holder, his (her) estate receives the 
remaining amount in the risky asset account. Typically, these contracts have ratchet provisions (step-ups), 
that periodically increase the benefit base if the risky asset investment has increased to a value larger than 
the guarantee account value, and roll up provisions, that periodically increase the benefit base according to 
a deterministic function. In addition, the benefit base may also be increased if the contract holder does not 
withdraw in a given year (bonus). Finally, the contract holder may withdraw more than the contractually 
specified amount, including complete surrender of the contract, upon payment of a penalty. Gomplete 
surrender here means that the contract holder withdraws the entire amount remaining in the investment 
account, and the contract terminates. In most cases, this penalty for full or partial surrender declines to 
zero after five to seven years. 

The hedging costs for this guarantee are offset by deducting a proportional fee from the risky asset 
account. From an insurance point of view, these products are treated as financial ones: the products are 
hedged as if they were pure financial products, and the mortality risk is hedged using the law of large 
numbers. Therefore, it is very important for insurance companies to be able to price quickly these products. 
Moreover these products have long maturities that could last almost 60 years. The Black-Scholes model, 
with constant interest rate and volatility seems to be unsuitable for those products: that’s why we present 
our pricing methods in two frameworks, modeling stochastic volatility (Heston model [TT]) and stochastic 
interest rate (Hull-White model [TT]) . 

There have been several recent articles on pricing GLWBs. In particular, we would remember the Forsyth 
and Vetzal’s work |3]: they used a PDF approach in a multi regimes model to price GLWBs contracts. This 
approach proved to be very fast and accurate, and we used it as a reference for our work. Goncerning 
the use of stochastic volatility, Kling et al. m used a Monte Garlo approach to price products. We have 
made reference also to Bacinello et al. [3: variable annuities (including GLWBs) are priced using a Monte 
Garlo approach. The policy holder (hereinafter, we will abbreviate it with PH) behavior is assumed to be 
semi-Static, i.e. the holder withdraws at the contract rate or surrenders the contract. 

In this paper, we price GLWBs guarantees, and we find the no-arbitrage fee, in the Heston model and 
the Black-Scholes with stochastic interest rate model {BS HW model). First, we treat a static withdrawal 
strategy: the PH withdraws at the contract rate. Then, taking the point of view of the worst case for the 
hedger, we price the guarantees assuming that the contract holder follows an optimal withdrawal strategy. We 
also used these methods to calculate the Greeks for hedging and risk management. Moreover we performed a 
mortality shock useful in risk management framework. For this purpose we present four numerical methods: 
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a hybrid tree-finite difference method and a hybrid Monte Carlo method (both introduced by Briani et al. 
[3]) an ADI finite difference scheme (Haentjens and Hout [13]), and a Standard Monte Carlo method with 
Longstaff-Schwartz least squares regression (Longstaff and Schwartz |16|h 

We use the term no-arhitrage fee in the sense that this is the fee which is required to maintain a replicating 
portfolio. A description of the replicating portfolio for these types of guarantees is given in Chen et al. |5] 
and Belanger et al. [3|. 

The main results of this paper are the following ones: 

• We formulate the determination of the no-arbitrage fee (i.e. the cost of maintaining a replicating 
hedging portfolio) in the Heston model and in the BS HW model using different pricing methods; 

• We present the effects of stochastic volatility and stochastic interest rate on pricing and Greeks calcu¬ 
lation, and the sensitivity of the GLWB fee to various modeling parameters; 

• We use different numerical methods to price the GLWB contract; 

• We present numerical examples which show the convergence of these methods. 

The paper is organized as follows: in Section 2, we describe the main features of the contract such as 
mortality, withdrawals, and ratchets. In Section 3, we provide a brief review of the stochastic models used 
afterward. In Section 4, we present the numerical methods, and how to implement them to solve the GLWB 
contract pricing problem. In Section 5 we perform tests in order to show their behavior and we study the 
sensitivity of the no-arbitrage fee to economic, contractual and longevity assumptions. Finally, in Section 6, 
we present the conclusions. 


2 The GLWB Contract 

In the following, we will refer to the contract described in the paper of Forsyth |3|, with some variations 
useful to compare our results with other works. We make a brief summary of the main features of the 
contract. 

2.1 Mortality 

We price the products in a risk-neutral measure, therefore in the following we assume that mortality risk 
is diversifiable (Milevsky and Salisbury, (TT]). When this assumption is not justified, then the risk-neutral 
value of the contract can be adjusted using an actuarial premium principle (Gaillardetz and Lakhmiri, |1()|L 
Hereinafter, the time variable will be denoted by the letter t, and we assume that the contract starts at 
t = 0. 

First we suppose that no PH can live longer than a given age. This age will be denoted by r (usually 
T = 122). The age of the PH at the beginning will be denoted by oq (usually oq = 65). So, the maturity of 
the contract is T = t — gq (usually T = 57): when the time variable t reaches T all PHs are died, and the 
contract is worth zero. 

The effects of the mortality on the contract are described using two functions: 

• Af : [0, T] —>■ K is the probability density that describes the random variable M associated to the death 
year of the PH. The fraction of the original owners who die in [t, t -I- dt] is equal to Af (t) dt. 

• 77. : [0, T] —K is the fraction of the original owners who are still alive at time t 

77(7) = 1— f M{s)ds. 

Jo 

We remark that 77 (0) = 1 and 77 (T) = 0. For seek of simplicity, we assume Af to be constant between 
contract’s anniversaries: if 7 S [fc, fc -I- 1[, k gN then Al (7) = A4 (k). 
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2.2 Contract State Parameters 


At time t = 0 the policy holder pays with lump sum the gross premium GP, to the insurance company. This 
may be rednced by some initial fees, giving a net premium P. The premium P is invested in a fnnd, whose 
price is denoted by the variable St- The state parameters of the contract are: 

• Account value: At, Aq = P. 

• Base benefit: Bf, Bq = GP. 

Both these two variables are initially set equal to the gross premium or to the premium. 

We suppose that the acquisition charges are equal to GP — P aren’t nsed for hedging purposes, but only 
to cover entry costs for management control. We suppose that there is a set of discrete times ti, which we 
term event times. At these times, withdrawals, ratchets, and bonnses may occur. Normally, event times are 
annually or quarterly. We first consider the evolution of the value of the guarantee excluding these event 
times ti. 

The value of the contract at time t is denoted by V {At, Bt,t). 

2.3 Evolution of the Contract between Event Times. 

Let t G ]ti,ti+i[ C [0,r]. As we said before, St denotes the underlying fund driving the account value. The 
dynamics of St will be described in the next Section. The acconnt value At follows the same dynamics of St 
with the exception of the fact that some fees may be snbtracted continuonsly: 

dAt = ^dSt - atotAtdt. (2.1) 

We snppose that total animal fees are charged to the policy holder and withdrawn continnonsly from the 
investment account At. These fees include the mutual fund management fees am and the fee charged to 
fund the guarantee (also known as the rider) ag, so that 


The only portion used by the insurance company to hedge the contract is that coming from a^: the other 
fees has to be considered as a ontgoing money flow as PH’s withdrawals are. 

Continuonsly withdrawn fees are typical of the contract described by Forsyth. Fees may also be withdrawn 
at the end of each policy year ti'. this is what Kling et al. do in m- In this second case 

dAt = ^dSt. ( 2 . 2 ) 

at 

When the PH dies, a death benefit, nsually equal to At, is payed out to the heirs of the PH. According to 
the contract’s formulation, this death benefit may be payed immediately or at the npcoming event time. If 
it is payed immediately, the contract stops immediately and the acconnt value and the benefit base becomes 
eqnal to zero; otherwise the contract goes on up to the next event time as if nothing happened. 

2.4 Event Times 

An event time is a seqnence of operations nnder the contract, which occnr at fixed dates, nsnally at each 
anniversary of the signing of the contract. The times these events take place are denoted hy ti — At ■ i and 
usnally At = 1. Let’s define I = T/At, then, i runs in {0,... /}. 

When an event time occurs, we assume that the following events happen in this order: 

1. Withdrawal of the fees by the insurance company (if it is not time continuons); 

2. If the PH died, payment of the death benefits; 
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3. If the PH is still alive, he (she) is entitled to withdraw a certain amount of money; 

4. If provided by the contract, a ratchet may increase the benefit base Bf. 

We denote with (H)(, Hj”, t^) the state variables just before an event time that occurs at time ti and with 
,ti) the state variables just after the update due to the i-th point of the previous numbered list. 


2.4.1 Fees 

Fees may be withdrawn continuously by the account value, as supposed in Forsyth and Vetzal injH]. In this 
case, between two event times, the account value changes as prescribed by ( 2 . 11 , and nothing special happens 
at an event time: 

Otherwise, fees may be withdrawn at the end of the period, as supposed in Kling et al. [15]. In this case, 
between two event times, the account value changes as prescribed by ( | 2 . 2 [ ), and at the event time, the account 
value becomes 

{Al+,Bl+,u) = . 

It is important to be able to deduce the management fees F™®” withdrawn by the account value because 
they are not used to hedge the contract and therefore they have to be considered as an outgoing money flow. 
If these fees are withdrawn continuously, we can calculate them observing that their dynamic between two 
event times is 


= a^Atdt + ndt. 


This ODE has the following solution 


rr<man 


= / e 


j: 


,, du 


Oirn-^sds. 


and can be used in a Monte Carlo approach. 

If the fees are withdrawn at the end of the period, we can calculate management fees as a fraction of the 
total fees withdrawn: 

+ Al (1 - , 


T^man _ rpman i 


^tot 


' p'tot _ jptot 


2.4.2 Death Benefit 

If the policy holder died at an instant t S his (her) heirs will obtain a death benefit, that is usually 

equal to the account value. If the contract provides that the death benefit is paid immediately, then the 
death benefit DBt is paid in t and is equal to At- Otherwise, if the DB is payed at the next event time, 
DBt^ = Al^ and the contract is concluded (after the DB payment it’s worthless): 

{A‘i+,B^+,U) = (0,0,F)- 

2.4.3 Withdrawal, Bonus, Surrender Event 

According to the contract, the policy holder, if still alive at event time ti, is entitled to withdraw a certain 
amount Wt, from his (her) police, also if the account value is equal to 0. This amount is given by 

Wu = GAt ■ B^+, 

where G is a constant defined by the contract. In a static framework, we suppose that the PH simply 
withdraws WAt^. Otherwise, in a optimization framework, he (she) may withdraw a fraction 7 ^ of the 
guaranteed withdrawn: 

Wt, = 7zGAf • Bl+. 


5 




• The case 7 ^ = 0 corresponds to no withdrawal. In this case, the contract may provide a bonus (&*, is 
specified by the contract): 


k) = {A^+, B^+ (1 + bt,), k) . 


• If 0 < 7 i < 1 the PH withdraws at a lower rate than the standard rate, and the new state variables are 

= (max (O, Al+ -Wu) , B^+,U) • 

• A third case is possible: the PH may want to withdraw more than the maximum admitted. In this 
case we suppose 7 ^ S ]1, 2], where the case 7 i = 2 corresponds to a total surrender. We define 


A' = max ( 0 ,A?+-GAt.H 2 +). 


The withdrawn amount is 


= GAt ■ B^+ + ( 7 , - 1) A' (1 -kJ. 

where Kt- € [0,1] is a penalty for withdrawal above the contract amount. The new state variables are 
(Al+,Bl+, t,) = (max (O, Al+ - GAt ■ B^+ - ( 7 , - 1) A'), (2 - 7 ,) B^+,u) 

= ((2-7,)A',(2-7,)H2+,t,). 


2.4.4 Ratchet 

If the contract species a ratchet (step-up) feature, then the value of the benefit base B is increased if the 
investment account has increased. The guarantee account B can never decrease, unless the contract is 
partially or fully surrendered: 


{At +, Bt +, k) = {Al+ , max {B^+ , A^+) , u) . 


Another feature that may be included in the contract is roll-up: for seek of simplicity we won’t treat this 
mechanism. 


2.5 Similarity Reduction 

An important property of GLWB contract is the fact that these contract behave good under scaling trans¬ 
formations. If V {A, B, t) denotes the price of a contract, it is possible to prove that for any scalar 7 > 0 


77 V {A,B,t) = V {r]A,r]B,t) ■ 


(2.3) 


Then, we just have to treat the case B = B for a fixed B (for example B = P), and then, choosing 
r] = ^js, we can obtain 

V(A,H,t) = |v (^A,B,^ , 

which means that we can solve the pricing problem only for a single representative value of B. This 


effectively reduces the problem dimension. The similarity reduction (2.3) was also exploited from Shah et 


Bertsimas in m- We can observe how the reduction similarity works both in the case of a contract that does 
not contain mechanisms for increasing the base benefits (ratchet), both for contracts with these properties. 


3 The Stochastic Models of the Fund S 

To understand the different impacts of stochastic volatility and stochastic interest rate over such a long 
maturity contract, we price the GLWB VA according to two models: the Heston model, which provides 
stochastic volatility, and the the Black-Scholes Hull-White model, which provide stochastic interest rate. As 
we said before, the process S represents the underlying fund driving the product’s account value At- 
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3.1 The Heston Model 


The Heston model m is one of the most known and used models in finance to describe the evolution of the 
volatility of an underlying asset and the underlying asset itself. In order to fix the notation, we report its 
dynamics: 

idSt=rStdt + VVtStdZt^ So = So, 

\dVt = k{9- Vt) dt + uj^tdZY Vo = Vo, 

where Z^ and Z^ are Brownian motions, and d {Z^, ZY) = pdt. 

3.2 The Black-Scholes Hull-White Model 

The Hull-White model |T1] is one of historically most important interest rate models, which is nowadays 
often used for risk-management purposes. The important advantage of the HW model is the existence of 
the closed formulas for prices of bonds, caplets and swaptions. In order to fix the notation, we report the 
dynamics of BS HW model: 

f dSt = nStdt + aStdZf = So, 

Idrt = fc {6t - rt)dt + ujdZ^ tq = fo, 

where Z^ and Z^ are Brownian motions, and d(Zf, ZY) = pdt. 

The process r is a generalized Ornstein-Uhlenbeck (hereafter OU) process: here 9t is not constant but it 
is a deterministic function which is completely determined by the market values of the zero-coupon bonds by 
calibration (see Brigo and Mercurio [7]): in this case the theoretical price of ZCB match exactly the market 
prices. 

Let (0, T) denote the market price of the zero bond at time 0 for the maturity T. The market 
instantaneous forward interest rate is then defined by 

It is well known that the short rate process r can be written as 

n = ujXt + Sit), 


where X is a stochastic process given by 


dXt = -kXtdt P dZl, Xo = 0, 


and S (t) is a function 

P it) = (0, ^) +(-^0 )^ ■ 

Then, the BS HW model is described by 


r dSt = nStdt + aStdZY So = So, 
< dXt = -kXtdt + dZl Xo = 0, 
[rt =UjXt + Sit). 


(3.2) 


A particular case is called flat curve. In this case, we assume it,T) = e and (0,T) = tq. 

Then 

2 

Sit) =ro + ^il- exp i-kt)) , 

and 


6*4 = -To + (1 - exp i-2kt)). 
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4 Numerical Methods of Pricing 

In this Section we describe the four pricing methods: an Hybrid Monte Carlo method, a standard Monte 
Carlo method, a Hybrid PDE method, and an ADI PDE method. 

We remember that our aim is to find the fair value for ag\ it’s the value that makes the initial value of 
the policy equal to the initial gross premium. To achieve this target, we price the police (with one of the 
following procedures) and then we use the secant method to approach the correct value for Ug. Therefore, 
the main goal is to be able to find the initial value for a given value of V (Ag, Bg, 0) (ag). 

We remark that we want to calculate the value of the police from the point of view of the insurance 
company: the management fees are treated as a outgoing cash flows, and if we assume that the policy holder 
follows a withdrawal strategy, we consider the worst one for the insurance company. 


4.1 The Hybrid Monte Carlo Method 

The value of a GLWB police can be calculated through a Monte Carlo set of simulations. This procedure 
is based in two steps: generation of a scenario (a sampling of the underlying values along the life of the 
product), and projection of the product into the scenario. According to the way we obtain the scenarios, we 
distinguish two Monte Carlo models: hybrid MC (HMC) and standard MC (SMC). 

The hybrid MC method was introduced in Briani et al. [S]. It is a simple and efficient way to produce MC 
scenarios for different models. This method is called “hybrid” because it combines trees and MC methods. 
First, a simple tree needs to be built: this can be done according to Appolloni et al. and also m , or as 
we are going to show in Em Then, using a vector of Bernoulli random variables, we move from the root 
through the tree, describing the scenario for the volatility or the interest rate. The values of the underlying 
at each time step can be easily obtained using an Euler scheme. 


4.1.1 Trees 


The trees for the Heston model and BS HW model can be obtained from Appolloni et al. [5] or Nelson 
and Ramaswamy m- In this case, the trees are simple binary trees: the node values, and the transition 
probabilities are set in order to match an approximation of the first two moments of the processe. This 
kind of tree perform well on short maturity, but the approximation errors accumulate on long maturities. 
Because of this error that accumulates, the convergence of the algorithm proved to be slow. Therefore, it 
was necessary to rethink the trees: the main aim was to set up trees which matched exactly some moments 


of the processes to be diffused. Here we present two trees (see Figure 4.1), one for stochastic volatility and 
one for stochastic interest rate. They are simple quadrinomial trees, and they are built to match the first 3 
moments of the stochastic processes. 

We suppose to fix a number N > 0, and we define h = t'/n. 


The General Case Let Z he & Brownian motion, and let G be a Gaussian stationary process, following 


dGt = a {Gt) dt + bdZt, 

with variance that depends only by the time lapse, i.e. Gt+s — GsjJ^s Af (/i {t, Gs), <J^ (t)). We show how 
to build a simple quadrinomial tree that can match the first three moments. 

We define a quadrinomial tree. Let’s fix a maturity T, and the number of steps N. Each node will be 
denoted by G(^n,j) where n runs from 0 to N, and j from 0 to 3n. Let h = t/n. The value of each node is 

= Go + (j - 1.5n) \/cr2 {h). 

We remember the first three moments of the process G: 


Mi=E[Gt+^-Gt|J-t] 


M (^j Gt) , M2 


= E 
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{Gt+h 


Gt)^ \Tt 


{h, Gt) + (h ), 









1 ^ 


^{n.j 

Node (n.j) 


37A ll — cxp{- 2 kh) 


Node (3,8) 



Time t 


Figure 4.1: The trees for Heston and Hull-White models. 


= E 


{Gt+h 


Gtf At 


= A {h,Gt) + 3fi{h,Gt)a^ (h). 


Let’s fix a node G(^n,j)- To be brief, fj, will denote fj, (^h,G(^n,j)) and a will denote y/(h). We suppose 
that the expected value ^ falls between the values of the nodes at time {n + 1) h . This hypothesis can be 
obtained assuming that the time step h is small enough. 

We define 


JA (nj) = ceil 


Gq — ^ 


-f 1.5(n-f 1) 


i.e. the first node in the next time step level whose value is bigger than the mean of the process. This can 
be seen in Figure 4.1 (both sides): the arrow points out to the expected value of the process, and Ja {n,j) 
is marked on the Figure. Let 


JB {n,j) = Ja {n,j) - 1, jc (nj) = Ja {n,j) + 1, Jb {n,j) = jo {n,j) - 2. 


To be brief we will only write Ja, Jb, jc^ jo-, and Ga will be Ga = and the same for the other 

letters: this is clear in Figure |4T| on the right side. 

We can now define a Markovian discrete time process n = 0,..., N with Go = G(o,o) and we suppose 
that if Gn = G(n,j)^ then it can move to Ga, Gs, Gc, Gd, according to the following probabilities 


PA 


P 


G„+i = Ga|G„ = G(, 


■d) 


(Ga - p) {{Ga - cr - pf + A) 

2(7^ 


Pb = P 


Gn+l — Gb|G„ — G(^n,j) 


(/i — Ga + 0") ((Ga — + cr^) 




PC = P 


Gn+l — Gc\Gn — Gfjij) 


(p — Ga + 0 ") ((Ga — cr — p)^ + 2tT^) 


6a^ 


PD = P 


Gn+l — Gd|G„ — G(^n,j) 


2a^(GA — p) + (Ga — p)^ 

6a^ 


And since Ga — cr < p < Ga , we can easily show that these probabilities are well defined: all in [0,1], 
their sum is equal to 1 , and the first three moments of the variable G„+i|G„ = G(„j) are equals to the first 
three moments of the variable Gt+h\Gt = Gf^nj)- 

Now, we approximate the process G by a discrete process G that is constant in each time lapse, and is 
defined as Gt = G[t/jvj • The weak convergence of this tree can be proved as in Nelson and Ramaswamy |18| . 
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o 


Figure 4.2: The possible combinations used to get positive probabilities in the Heston model tree. The red 
points correspond to the nodes used. 


The Heston Model The Heston process (3.1 1 for volatility has no constant variance and isn’t Gaussian. 


We consider the process obtained by the square root: 

Ak (0 - 


d^/Vt = 


8 v^ 


-dt + —dZt- 


We approximate it with a Gaussian process with variance ^dt. This approximation is helpful to define 
the grid: inspired by m, we define 

jn = max ( 0, floor ( 1.5n- 

V V wV/i/ 


and we set 


V(„j) = ^max ^0, y/Vo + (j + jn - l-5n) 


UJ 


y/h' 


for j = 0,..., 3n — The shift due to helps to reject the many node with value equal to zero: if 
jn > 0, then V(^n,o) = 0 and V(^n,i) > 0 . 

We fix now the value of n and j. The discrete process V can jump from a node to another, as in a 
Markovian chain. We show now how to find the possible upcoming nodes. 

The first three moments for the Heston process can be found in Alfonsi [T]: 

{h) = Ml = E [Vt+h\Vt =v]= [h ), 


M2 = IE 


(Vt+hY\Vt=v = Ml + (h) [eki; {h) /2 + ve->^^] , 


M 3 


E 


(Vt+hY ir* = V 


M 1 M 2 + (h) 


2^2g-2/ch ^ 



(3ne-'=^ + ekij (h)) 


Then we can proceed as in the general case. Anyway, the grid we’re using is based on an approximation: 
so the probabilities obtained solving the linear system may not be positive. 

If we get negative probability for a given node, we try another combination of nodes: the node A or C 
may be replaced by a node E define as the first node bigger than C, and the node B or D may be replaced 
with a node F, defined as the smallest before node D. This gives rise to 9 combinations to be tested. If the 
starting node is small and the node D verifies jo = jn we could not do this last change because there would 
be no F node. In this case we allow the node D to be replaced by the node E: see Figure [T^ 

If these attempts don’t give a positive result (negative probabilities), we give up trying to match the first 
three moments, and we are content to match an approximation of the first two as in |18| . thus ensuring the 
weak convergence. In this case, we only use the nodes A, B, C, D: we define 

M V^n+l,jB 1 

PAB = ^ PBA = 1 - PAB, 

^n+ljA ^n+l,jB 
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PDC = 1 - PCD 


PCD = 


Gn 


+ 1JC 


-G 


ra+l J'd 


and 


PA = -^PAB, Pb = -^Pba, pc = -^PCD, Pd = '^Pdc- 


It is possible to show that the first moment of this variable is equal to Mi, and as h —>■ 0 the second 
moment approaches to M 2 , ensuring the convergence, as proved in |18| . 

In all our numerical tests, this last option (matching only two moments) has never been necessary: 
changing the nodes, all moments were matched with positive probabilities. 


The Hull-White Model The process X in (3.2 1 is Gaussian. As shown in Ostrovski |20] the variables 
Xt and Xydy are bivariate normal distributed conditionally on X^ with well known mean and variance. 
We define 


j) ~ [J 2 


1 — exp (—2kh) 


2k 


, n = 0,..., N and j — 0,..., 3n. 


Let’s fix a node A(„ j). We define 

H = exp {—kh), K = 


1 — exp {—2kh) 


, Ml = X(^n,j)H, 


jA = ceil 

The transition probabilities are given by 


Ml 3 (n + 1) 
~K 2 


, Xa — X(n+l,jA)- 


PA = + (K + Mi- XAf) , pb = (x^ + (Ml - XAf) , 

PC = (2K^ + (K + Mi- XAf) , PD = (2K^ + (Mi - XAf) . 

4.1.2 Scenario Generation 

The generations of the volatility process and of the interest rate process behave in a similar way: we start 
from the node (0,0) of the tree and according to a discrete random variable and to the node probabilities, 
we move to the next node and so on. Let D be a discrete random variable that can assume value A, B, C, D 
with probabilities pa,Pb,PCiPd- sampling such a variable at each node, we get the values of the process at 
each time step. 

We distinguish two cases for the two models. 

The Heston Model We approximate the couple (S'*, 1/*) in [0, T] by a discrete process {SkAt,VkAt) t/ad 
with (5o,h()) = (5'o)hb)- For each scenario, we generate the volatility. 

Let iV ~ Af (0,1) and B ^ B (0.5). We deduce the value of St+At by 


St+At — 



(r - Pkd)At+(Pk - At+^(Ft+At - Vf + ^(l-p2)AtVtN 

(r - §k0)At+(§k - l) h+%^+^j At+§(Vt+At - Vt) + V(l - P^) ^tVt+AtX 


ifB = 0, 
ifB = 1. 


According to (13, we use the normal variable N to generate the Gaussian increment of S, and the 
Bernoulli variable B to split the operator associated to the Heston process. 

This scheme (without splitting) appears in Briani et al. |5] and the splitting method appears in Alfonsi 

ra¬ 


il 






















The Black-Scholes Hull-White Model We approximate the couple {St,Xt) in [0,T] by a discrete 
process (SkAt, XkAt)f.^Q t/av “ (*^ 0 ’ deduce the interest rate by ft = ujXt + (3 (t). 

Let N Af {0, 1). We deduce the value of St+At by 


St+At = St exp 


f rtAt 


--T'" 


(j^Xt+At + Xt {kAt — 1)) p + VAtpN^ 


4.1.3 Projection 

Once we have generated the scenarios, we project the police into it: it means we calculate the initial value of 
the contract as the sum of discounted cash flows. This calculation depends on whether we take an optimized 
strategy or not. Let V be the value of a police at time t, having account value equal to A and base 

benefit equal to B. From now on, we fix a specific scenario. Let V {A,B,t)i>G the value of a police in that 
scenario, at time t, having account value equal to A and base benefit equal to B. 


Constant Withdrawal In this case the strategy of the PH is fixed: in each event time 7 ^ = 1 (for 
completeness we continue to write 7 i). A simple way to calculate the value of the police is calculating 
forward the cash flows, conditioning on the death time. As in Holz et al. m, we have: 

V{Ao,Bo,0)=J2m iu) (^e- + e" ) . 


2 = 0 


\k^0 


Anyway, we developed another approach, useful for the optimal withdrawal case. First we calculate the 
values for all ti neglecting the effect of mortality (equivalently, assuming that the PH die at 

the end), with a forward approach: 


Ay = max 


Bt = 


( 0 , - 7qGAtHy^) 

max j, Ay ^ if ratchet. 


B. 


otherwise. 


Then, we proceed backwards, calculating the value of the contract for each time ti just before the withdrawal. 
The value of the contract at time ti can be written as the discounted value at time plus the discounted 
value of the cash flows relating the period [ti~^,ti+i] ■ The final condition on the value of the contract is 

p(Ay,Hy,T) = 0, 

because all PHs are death and all benefits have been paid. Then 


p(Ay,Hy,y = e 




V 


1 j ^y 1 ! ^i+i) + (tz+i) kFti+i 


+ DB + MF, 


where DB and MF stand for the discounted value in ti of the death benefit and management fees paid in 
[ti^,ti+i]- We distinguish four cases depending on how the management fees and the death benefit are 
payed. 


CASE 1: DB payed at the end, Fees withdrawn at the end 

Ai+ njrp _ -rp UA 


DB = M (ti) e~ •'*. MF = 7^ (L) e 

CASE 2: DB payed at the end, Fees withdrawn continuously 


rgds ,, 4 + 

St, 


1 ( 1 -, 


t\ Otm 
/ atot ‘ 


= At (L)e--^*yG,d.^4+_ 


^C-+l ^-atotAf 
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CASE 3: DB payed immediately, Fees withdrawn at the end 


a 4+ rti+i 

DB ^ M{ti) 

Jt- 


dt. 


MF = M (U) f St (l- e”-''**. + ('l _ 

Oitot J t- ^ ^ ^ 


C^tot 


CASE 4: DB payed immediately, Fees withdrawn continuously 

4'^+ r^-i + l r-t 

Jti 


Proceeding in this way, it is possible to calculate P (Aq"'', Bg"'', O). The initial value of the police is 

P(Ag-,i3o-,0) =P(A4+,i?4+,0), 


Af+ Bi+i 


if the first withdrawal takes place at time t = ti, or 

y (A„- ,Bo,0)=V (A4+ , B ^+, 0) + AoGAtP 


if the first withdrawal takes place at time t = 0. Then we simply have to calculate the average of 
V ( Aq , Bq , O) among the simulated scenarios. 


Optimal Withdrawal In this case we suppose that at each event time U the policy holder can withdraw 
a fraction 7 ^ of the regular amount. To price in this case, we suppose that the PH chooses the value of 7 
that causes the worst hedging case for the insurance company. We denote V {A,B,t) the expected value at 
time t of a generic police whose state parameters are A, B : 

V{A,B,t)=E[V (A, HA)]. 

So, we suppose that the policy holder chooses 7 ^ such that 

7 * = argmax V (A^+, 5"^+, t) . 

7e[o.2] 

This expected value can be calculated with a Longstafi-Schwartz approach: 

1. Simulate N random scenarios and price the police into these scenarios using random values for 7 ^. 

2. For t = r to t = 0: 

(a) Approximate the function V {A,B,t) using the least squares projection into a space of functions 
(usually polynomials). 

(b) For each scenario find the optimal value of 7 * . 

(c) Recalculate the upcoming state variables from s = t to s = T assuming that the policy holder 
chooses the best value for 7 . 

3. Calculate the average of the initial value V (Ag, Bq, 0) for all the scenarios. 

The approximation of the function V (A, H, t) can be improved by the reduction property. 
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4.2 Standard Monte Carlo Method 


The Monte Carlo method is very similar to the hybrid Monte Carlo one. The only different thing, is the way 
we produce the random scenarios. The projection phase is the same as hybrid Monte Carlo’s one. 

4.2.1 Scenario generation 

We distinguish two cases for the two models. 

The Heston Model The generation of the scenarios (underlying and volatility) in this case has been done 
using a third order schemes described in Alfonsi [1]. 

The Black-Scholes Hull-White Model The generation of the scenarios (underlying and interest rate) 
in this case has been done using an exact schemes described in Ostrovski |^, with a few changes in order 
to incorporate the correlation between underlying and interest rate. 


4.3 PDE Hybrid Method 

The Hybrid PDE approach is different from the previous ones. In fact it’s a PDE pricing method and it’s 
based on Briani et al. 111,0 both for Heston and Hull-White case. Using a tree to diffuse volatility or 
interest rate, we freeze these values between two tree-levels and we solve four PDE for each tree node; then 
we mix the values given by the four PDEs according to the transition probabilities of each node. 

We can resume the pricing methods in three features: model, algorithm structure and pricing. 


4.3.1 The Heston Model 

Starting from the model 


dSt 

dVt 


we define the process 


rStdt + VVtSt {pdZY + pdZf) Vq = Vq, ^ q 

k{e-Vt)dt + ujVVtdzY So = So, \ 


Et = In (A) - ^Vt, Eo = In (Ho) - 

LU UJ 

At = exp (^Et + ^Vtj . 


Then 



if fees are taken continuously, otherwise 



^k{9-Vt)- ato?j dt + ^{1-p^)VtdZf, 


Vt)^ dt + py/(l - p2) VtdZY- 


4.3.2 The Black-Scholes Hull-White Model 

Starting from the model 

(dSt = nStdt + aSt {pdZl + pdZf) So = So, 

< dXt = -kXtdt + dZl Xo = 0, d {Zf, Zl) = 0, 

[rt = ioXt + P{t), 


(4.1) 
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we define the process 


Ut = In {At) - paXt, Uq = In (Aq) , 
At = exp {Ut + pcrXt). 


Then ^ 

dUt ^ - Y crpkXt - atot^ dt + a^/l^^dZ^, 

if fees are taken continuously, otherwise 


dUt = 



apkX^ dt + a\/l 


p^dZ't 


(4.2) 


4.3.3 Algorithm structure 

The structures for this algorithm consist in a tree and a PDE solver. As described in Briani et al. mM , we 
use a tree to diffuse the volatility (or the interest rate) along the life of the product, and we solve backward 
a ID PDE freezing at each node of the tree the volatility (or the interest rate). The tree is built according 
to Section 4.1.1 (quadrinomial tree, matching the first three moments of the process), and the PDE is solved 


with a finite difference approach. We have to solve the PDE between event time, and at the event time we 
apply the changes to the states to reproduce the effects of the events. 


4.3.4 Pricing 

The PDE we have to solve is essentially the same as in Forsyth and Vetzal [S]. We distinguish four cases 
as we did in Monte Carlo case. We denote with V {A, B, t) the value of a contract at time t, whose account 
value is worth A and whose base benefit is worth B . Consequently, we define 


and 


(E, B,t)=V (exp (e +^Vt), B, t) , 


{U,B,t)=V (exp {U + paXt) B,t). 


The variables f, X and V will denote the frozen values of rt, Xt and V*. We solve the transformed PDE 
between two event times for each node of the tree four times: one for each of the possible next nodes, using 
the initial data corresponding to these nodes. To reduce the run time, we do this only for most relevant 
nodes: this cutting technique dramatically reduced calculation times without compromising the quality of 


results. Then, using the inverse transformations (4.11 and (4.21, we apply the event times’ actions. In the 
next few paragraphs, we are going to write 2 PDFs: one for the Heston model, and one for the BS HW 
model. 


CASE 1: DB payed at the end, Fees withdrawn at the end 

The terminal condition is 


V{A,B,T) 


TZ{T- At) A 


(^1 - (1 - 



The associated PDEs are 


vf" + + (^ - I - - ^)) = 0 ’ 

—2 2 / 2 \ 

IjHW , P C . I - ^ , „uv\ yjHW n 

Vt H-V T ^ o-pkX j Vu -ry = 0 . 

For ti = T — 1 to ti = 0 we have to: 


(He 1) 
(HW 1) 
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1. Solve the PDE backward from ti+i to U. 

2. Calculate the value of V from the value of or 

3. In case of ratchet V [A, B, tf'*') = V {A, max {A, B ), . 

4. Withdrawal: 

(a) if = 0 : 

V{A,B,tl+)=V{A,B{l + bt,),i*)- 

(b) if e [ 0 , 1 ] : 

V {A, B, tl+) = V (max (0, A - 7 *. GAtB ), B, t®+) + 7^ (b) 7 *, GAtB; 

(c) if 7 t. e ] 1 , 2 ] : 

V (A, B, t?+) = V (max (0, A - GAtB) (2 - jtJ , B (2 - jt ,), t?+) + 

+ Ti{ti) {GAtB + {yti — 1) max (0, A — GAtB) (1 — Kt J); 

5. Death benefit: V {A, B, t\+) = V {A, B, t^+) + (7^ (b_i) - 7^ (U)) A. 

6 . Fees: V (A, B, t") = V (^ 6 -“*°*^^*, B, ti+) + 7^ (b_i) (l - 

7. Calculate the value of or from the value of V . 

CASE 2: DB payed at the end, Fees withdrawn continuously 

The differences between this case and the case 1 are the following ones. The terminal condition is 

V{A,B,T) =n{T-At)A. 

The associated PDFs are 

Vf" + + (^ “ f “ 5^ + OLmB (t) exp (b* + gp) = 0 , (He 2) 

^Hw ^ - ato^ - fV"^ + a^TZ (t) exp (Ut + paX) = 0 . (HW 2) 

Point 6 (fees step) becomes 

V{A,B,t-)=V{A,B,tl^). 

CASE 3: DB payed immediately, Fees withdrawn at the end 

The differences between this case and the case 1 are the following ones. The terminal condition is 

V{A,B,T)=0. 

The associated PDFs are 

Vf" + ^Vf|+(^r-|-^fc(e-P))vf"-rV^"+Af (b)exp(Bt + ^p) (^ 1 -= 0, (He 3) 

V!^^ + ^V^^+(r-^+apkX^ (b)exp(C*+paX)(^l-(l-e-“‘°‘(‘-‘-)) ^) = 0 • (HW 3) 

Point 5 (death benefit step) and 6 (fees step) become: 

• Death benefit: V (A,B,ti^) = V (A, B,t^^)- 

• Fees: V {A, B, t") = V (Ae-“‘°*^*, B, tl+) + 7^ (b) (l - . 
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CASE 4: DB payed immediately, Fees withdrawn continuously 

The differences between this case and the case 1 are the following ones. The terminal condition is 

V{A,B,T) = 0. 

The associated PDEs are 

Vf' + ^Vee + (^ - f + exp (^Et + (amTZ (t) + M {U)) = 0 , (He 4) 

yHw ^ - atot) + exp {Ut + paX) (a™7^ (t) + M {U)) = 0 . (HW 4) 

Point 5 (death benefit step) and 6 (fees step) become 

V (A, B, t-)=V {A, B, tl+) = V {A, B, . 


This concludes the static withdrawal case. In the optimal withdrawal case, we suppose the PH to change 
the value of 7 i used in step n. 4 (withdrawal step). He (she) will choose the value of 7 ^ G [0,2] in order 
to maximizes the value of V (A,B,ti~^)- This maximization can be done using a grid of values for ji and 
choosing at each time the best value. 

4.4 PDE ADI Method 

Consider the asset price process given by the system of stochastic differential equations described in Section 
We describe the ADI method only in the case 2, but the other cases can be easily adapted. Moreover, we 
have chosen to not use the transformed PDE described in Section |4.3.4[ but the classical version of PDEs 
for the Black-Scholes, Heston and Black-Scholes Hull-White model. The associated PDEs are 

vf" -t + ir- atot) AVa" + PujAVVav + k {6 - V) Vy" - rV"" -f a,^TZ {t) A = 0 (He 2b) 


Vf'^ + 


VlT -t + {r- atot) AvS'^+ pujAaVA^ + k{9t- r) Vf- rV"'^ -f Q^7^ (t) A = 0 (HW 2b) 


Because of the long maturity, solving a two-dimensional PDE is a very costly and slow method. The idea 
is to use splitting schemes of ADI (alternating directional implicit) type. In this paper, we only present the 
Douglas scheme, but various scheme are available in the literature. In order to solve the PDE, we should 
address many numerical difficulties. The first one is the mesh and we have chosen to use the meshes described 
in m with the parameters 


Aieft = C.BiSo Aright = l-2S'o Amax = 1005'o and di = 5'o/20, 


for the mesh of variable A, 


I^max — 10.^0 5 ^ — ^0 and C?2 — .^max/400 

for the mesh of variable r in the Black-Scholes Hull-White model, and 

Kiax = MIN{MAXimVo, 1), 5) and = Enax/500. 

for the mesh of variable V in the Heston model. The second difficulty is the choice of the splitting scheme. 
We have chosen the Douglas scheme with parameter 0 = 1/2 because it is the easiest to implement, but of 
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course some higher order schemes (in time) would be more optimal. The last difficulty, but not the least, 
is the choice of boundary conditions. Since there is no closed form solutions for the GLWB product, it is 
difficult to make the right choice for the boundary conditions. Moreover the boundary conditions have a big 
impact on the solution, because of the long maturity. The choice of homogeneous Neumann conditions is 
usually done because it simplifies the system to solve (exactly it simplifies the finite difference scheme at the 
boundary). In the context of GLWB, the boundary conditions for the Black-Scholes Hull-White model will 
be given by: 


gyHW 

(A r,t) = 0, \i A = 0 or A = Hmax, 

q^^hw 

{A, r, t) = 0, if r = ±i?max, 

on the mesh [0, Hmax] x [~.Rmax, .Rmax], and the boundary condition for the Heston model will be given by: 

{A, V,t)=0, if H = 0 or H = H^ax, 

-^{A,VA)=0, ify = K.ax, 
ov 

on the mesh [ 0 , Amax] x [ 0 : Knax]? with no condition at = 0 since it is an outflow boundary. 

5 Numerical Results 

In this Section we compare the numerical methods used in Section]^ Hybrid Monte Garlo {HMC), Standard 
Monte Garlo {SMC), Hybrid PDE {HPDE), and ADI PDE (APDE). In particular we compare pricing and 
Greeks computation in Static Case [5T] and Dynamic Case [5?^ 

We chose the parameters of the methods according to 4 configurations (A, B, C, D), with an increasing 
number of steps, and so that the calculation time for the various methods in each configuration were close. 
There 4 configurations are in Table with the notation (time steps per year ; space steps; vol steps) for the 
ADI PDE method, (time steps per year ; space steps ) for the Hybrid PDE method approaches, and (time 
steps per year ; number of simulations) for the MG’s one. In Monte Garlo for dynamic case, we also add the 
degree of the approximating polynomial. These values had been chosen to achieve approximately these run 
times: (A) 30 s, {B) 120 s, (C) 480 s, {D) 1900 s. To reduce the run time we do the secant iterations using 
an increasing number of time steps for all the methods: the values in Table [T] are those used for the last 3 
iterations. 

We use the standard MG both as a pricing method, both as a benchmark (BM). About the benchmark, in 
the static case we used 10^ independent runs. In the dynamic case we used 10® independent runs, arranged 
in 10 sub runs; in each sub runs the expected value has been approximated by a 6 order polynomial. At 
each event time, the PH can chose between 7 = 0, 7 = 1 and 7 = 2 . 

The search for the fair ag value has been driven by the secant method. The initial values for this method 
were = 0 bp and Og = 200 bp. 

To achieve Delta calculation in Monte Garlo methods we used a l%o shock in static case and 1% in 
dynamic case. 

We used the DAV 2004R mortality Table, 65 year old German male (see [3] for the Table). It contains 
the probabilities that a person aged t will die within the next year. It’s easy to get the function Ai from 
these probabilities. 

5.1 Static Case 

In the static case we suppose the PH to withdrawal exactly at the guaranteed rate: 7 ^ = 1. 
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BS HW Static 



Heston Static 



HMC 

SMC 

HPDE 

APDE 

HMC 

SMC 

HPDE 

APDE 

A 

5x1.3-10® 

1x2.7-10® 

30x400 

18x180x36 

5x8.6-10-’* 

5x7.4-10"* 

35x400 

26x260x13 

B 

10x2.3-10® 

1x9.8-10® 

60x600 

27x270x54 

10x1.6-10® 

10x1.4-10® 

70x600 

40 X 400 X 20 

C 

20x5.4-10® 

1x4.9-10® 

100x1000 

40x400x80 

20x3.8-10® 

20x3.5-10® 

100x1000 

64 X 640 X 32 

D 

40x1.0-10® 

1x2.0-10^ 

200x2000 

62x620x124 

40x7.3-10® 

40x7.5-10® 

200x2000 

104x1040x52 




BS HW Dynamic 



Heston Dynamic 



HMC 

SMC 

HPDE 

APDE 

HMC 

SMC 

HPDE 

APDE 

A 

5x3.3-10®x2 

5x3.2-10®x2 

30x400 

16x160x32 

5x3.2-10®x2 

5x3.2-10®x2 

35x400 

22x220x11 

B 

10xl.610"*x3 

5xl.6-10"*x3 

60x600 

24 X 240 X 48 

10xl.510"*x3 

10xl.510"*x3 

70x600 

36x360x18 

C 

20x5.2-10"*x4 

5x5.3-10"*x4 

100x1000 

38x380x76 

20x4.910"*x4 

20x4.&10"*x4 

100x1000 

60x600x30 

D 

40xl.410®x5 

5xl.6-10®x5 

200x2000 

60x600x120 

40xl.310®x5 

40xl.310®x5 

200x2000 

100x1000x50 


Table 1: Configuration parameters for the BS HW model and for the Heston model, static and dynamic. 


The static tests 1 and 2 are inspired by [^: in their article, Forsyth and Vetzal price a GLWB contract 
in a static framework, under the Black Scholes model with r — 0.04 and tr = 0.15. The contract parameters 
are reported in the Table the contract type corresponds to case 2 in Section [4T3] and [4T4| 


Initial age of PH 

65 

Gr. premium 

100 

DB payment 

next anniv. 

G 

0.05 

Initial fees 

0 

Ratchet 

Off/On (annual) 

Withdrawal rate 

First withdrawal 

1 per Y 

1 ^* anniv. 

Fees taken 

0 

cont.ly 

Strategy 

static (7 = 1 ) 


Table 2: The contract parameters for static tests (except Test 2B). 


They treated two cases: no ratchet, and annual ratchet. In the first case they get ag = 35.51 bp and 
in the second case ag = 64.92 bp. In Test 1 and Test 2 we introduce respectively stochastic interest rate 
and stochastic volatility to analyze the impact of these model developments on the fair guarantee fee. The 
parameters for interest rate and volatility models has been chosen to be plausible. 

To compare our results in the Heston model with Kling’s ones in m we performed test 2B. In this case, 
product parameters are reported in Table and correspond to case 1 in Section [4.1.3| and |4.3.4| 


Initial age of PH 

65 

Gr. premium 

100 1 

DB payment 

next anniv. 

G 

4.90%, 4.19% if ratchet 

Initial fees 

4% 

Ratchet 

Off/On (annual) 

Withdrawal rate 

1 per Y 


151 bp 

Strategy 

static ( 7 = 1 ) 

First withdrawal 

1 ®* anniv. 

Fees taken 

at the end 




Table 3: The contract parameters for Test 2B-Static. 


5.1.1 Test 1-Static: the Black-Scholes Hull-White Model 

In this test we want to price a product according to BS HW model. We use the same corresponding 
parameters as in test [S]. Model parameters are shown in the Table Results are available in Table 

All four methods behave well and in the configuration D, gave results consistent with the benchmark. 
HPDE proved to be the best: all configurations gave results consistent with the benchmark. Then APDE 
and SMC, and HMC gave good results too. SMC performed a little better than HMC: the first method 
simulates the underlying value and the interest rate exactly and so it is enough to simulate the values at 
each event time. HMC matches the first three moments of the BS HW r process, but doesn’t reproduce 
exactly its law: therefore it is right to increase the number of steps per year. So, for a given run time, we can 
simulate less scenarios in HMC than SMC: effectively, the confidence interval of HMC is larger than SMC’s 
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one. Moreover, SMC over performed the benchmark when using configuration D. Particularly, correlation 
between underlying and interest rate has a fundamental role, and it’s impact can be bigger than ratchet’s 
impact: for example, case no ratchet with p = 0.5 gave a higher price than case ratchet with p = —0.5 (111 
bp vs 84 bp). 

5.1.2 Test 2-Static: the Heston Model 

In this test we want to price a product according to the Heston model. Model parameters are shown in the 
Table [U Results are shown in Table 0 

In this Test, MC methods had more problems; PDE methods’ values are close to the benchmark, while 
MC method’s values were far, but compatibles with the benchmark (BM’s value is inside MCs’ confidence 
interval). Probably, in this case, the benchmark is not very accurate: this is due to the fact we used SMC to 
calculate it. If we compare the two MC approaches, in this case, they both use a third order approximation 
and than they become equivalent: HMC proved to be faster than SMC when using few time steps (we could 
exploit -1-16% simulations in configuration A), while SMC proved to be slightly faster in high time steps 
simulations, because of more time needed to build the volatility tree (—3% simulations in configuration D). 
HPDE showed to be very stable (case no ratchet, p = —0.5, Ug didn’t changed through configurations B-D), 
but APDE behaved well to (monotone convergence). In the Heston model, correlation has a less important 
role than in BS HW case: among the different values of p, the value of Ug changes less then 5 bp in no-ratchet 
case, and less than 1.5 bp in ratchet case. 

5.1.3 Test 2B-Static: the Heston Model 

In this test we want to obtain the results shown in m, where the contract are priced with MC techniques. 
The values given in m are 150 bp for both cases (no ratchet and ratchet). Model parameters are given in 
Table [8] Results are available in Table |9l 

In this Test, all methods gave the same results, but not the same results as in Kling et al. [15]. One 
possibility is that we have misinterpreted some of the contractual specifications in Kling’s paper, leading to 
some subtle differences in the contracts that we are considering as compared to theirs, and these discrepancies 
result in different fees. Another potential explanation is that a Monte Carlo method was used to determine 
the fee by Kling et al.; this may have introduced a significant error when calculating the fee unless a very 
large number of simulations was used. They didn’t report a confidence interval for their results, so it’s hard to 
understand the cause of the gap. Moreover we can observe that out two MC methods gave larger confidence 
intervals than Test 2-Static: probably, the parameters used for Test 2B-Static shape a harder pricing problem 
than previous test, and more simulations should be performed to obtain same quality results. Also in this 
case, HPDE proved to be the most stable method. 


5.1.4 Test 3-Static: Hedging 


To reduce financial risks, insurance companies have to hedge the sold VA: to accomplish this target they 
must calculate the greeks of products. 

In this test we want to show how the different methods can be used to calculate the main greeks. This 
can be done through finite differences for small shocks on the variable. In general, the PDE methods are 
ahead w.r.t. MC methods: the price is computed through finite differences and so the price under shock is 
already computed. For MC methods this is quite harder because the pricing has to be repeated changing 
the inputs. 

To start, we calculate the underlying greek delta, for the products of Test 1-Static and Test 2-Static. As 
in this case we don’t want to compute the fair fee ag, we fix it arbitrarily. We choose two values for each 
model: one for no ratchet case, and one for ratchet case. The values chosen are such as to cover the costs of 
the insurer regardless of the correlation, and may be plausible on a real case. Results are available in Table 
10 (all values in Table must be multiplied by 10“"^). 
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In this Test, we got very accurate results with all method. Anyway, HPDE proved to be the best: it is 
the more stable and accurate. We remark that despite fair fee changes a lot when changing the correlation 
parameter p, the value of Delta changes much less. Delta calculation proved to be harder in the Heston 
model case than in the BS HW model case. 

5.1.5 Test 4-Static: Risk Management 

Mortality and longevity risks are unhedgeable risks. Usually the Risk Management Team has to calculate the 
financial reserve taking into account these risks. Usually extreme scenarios are chosen and policies are priced 
according to them. In this test we analyzed how the different pricing methods behaved under mortality 
shocks: the mortality probabilities have been increased by 10% except the last one who’s equal to 1. To be 
brief, we simply report the fair fee for D case. Results are available in Table El 

In this Test, we got results similar to Test 1-Static and Test 2-Static, and mortality shocks didn’t affect 
the convergence quality of the four methods. We observe that mortality shocks reduce the value of ag (about 
minus 5 bp) and this means that an increase in mortality shouldn’t be a source of losses for the insurer. 
Consequently, insurers should pay attention to longevity risk. 
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Table 4: The model parameters about Test 1-Static. 




no ratchet 

annual ratchet 

P 


HMC 

SMC 

HPDE 

APDE 

BM 

HMC 

SMC 

HPDE 

APDE 

BM 

0.5 

A 

45.99 ± 1.06 

45.58 ± 0.73 

45.72 

47.35 

45.81 

84.85 ± 1.23 

84.58 ±0.85 

84.56 

88.20 

84.71 

B 

45.52 ± 0.79 

45.31 ± 0.38 

45.71 

46.09 


84.35 ± 0.91 

84.15 ± 0.44 

84.60 

85.66 


1 

c 

45.58 ± 0.52 

45.71 ± 0.17 

45.69 

45.99 

±0.12 

84.31 ± 0.60 

84.65 ± 0.20 

84.63 

85.36 

±0.14 


D 

45.85 ± 0.38 

45.71 ± 0.08 

45.72 

45.81 


84.68 ± 0.44 

84.63 ± 0.10 

84.64 

84.94 



A 

82.29 ± 1.39 

81.40 ± 0.95 

81.87 

82.91 

81.88 

157.77±1.65 

156.77±1.14 

156.36 

161.04 

157.09 


B 

82.43 ± 1.04 

81.53 ± 0.50 

81.92 

81.75 


157.68±1.23 

156.55±0.59 

156.46 

157.91 



c 

81.62 ± 0.68 

81.77± 0.22 

81.80 

81.89 

±0.16 

156.50±0.80 

157.05±0.27 

156.87 

157.70 

±0.19 


D 

81.99 ± 0.50 

81.83± 0.11 

81.79 

81.81 


157.16±0.59 

157.07±0.13 

156.96 

157.27 


d 

A 

111.75T1.76 

110.30±1.20 

111.14 

109.23 

111.05 

224.19±2.15 

222.14±1.48 

221.78 

227.14 

222.83 

B 

112.73±1.32 

110.85±0.63 

111.07 

108.93 


224.59±1.60 

222.26±0.77 

222.32 

223.44 


+ 

c 

110.89±0.86 

111.08±0.28 

111.05 

109.93 

±0.20 

222.18±1.05 

222.94±0.35 

222.52 

223.36 

±0.24 


D 

111.29T0.63 

111.11±0.14 

111.02 

110.42 


222.97±0.77 

222.94±0.17 

222.67 

222.96 




HMC 

SMC 

HPDE 

APDE 

A 

30 s 

30 s 

30 s 

28 s 

B 

119 s 

120 s 

128 s 

184 s 

C 

472 s 

478 s 

395 s 

461 s 

D 

1866 s 

1896 s 

1903 s 

1800 s 


2.5% 


2 . 0 % 


^ 1..5% 
o 
W 

1 . 0 % 


0.5% 


0 . 0 % 


Test 1-Static: Relative Errors (error/time) 


I 



0 500 1000 1500 2000 Time (s) 


- Cl 

— HMC 
— SMC 
— HPDE 
- • APDE 


Table 5: Test 1-Static. In the first Table, the fair fee for the Black-Scholes Hull-White model, with no ratchet 
or annual ratchet. In the second Table the run times for the no-ratchet case (p = —0.5). Finally, the plot of 
relative error (w.r.t. MC value) for the three methods in the case p = —0.5 with no ratchet. The parameters 
used for this test are available in Table |2] and in Table ID 
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Table 6 : The model parameters about Test 2-Static. 




no ratchet 

annual ratchet 

P 


HMC 

SMC 

HPDE 

APDE 

BM 

HMC 

SMC 

HPDE 

APDE 

BM 

0.5 

A 

36.77 ± 1.25 

36.17 ± 1.36 

37.00 

37.41 

37.16 

61.47 ± 1.23 

60.90 ± 1.35 

61.51 

62.30 

61.84 

B 

36.74 ± 0.92 

36.40 ± 0.99 

37.01 

37.26 


61.29 ±0.90 

60.85 ± 0.97 

61.59 

62.06 


1 

c 

36.79 ± 0.59 

36.94 ± 0.62 

37.01 

37.11 

±0.12 

61.36 ± 0.58 

61.56 ± 0.61 

61.63 

61.80 

±0.11 


D 

37.47 ± 0.43 

37.33 ± 0.42 

37.01 

37.06 


62.15 ± 0.42 

61.95 ± 0.42 

61.66 

61.77 



A 

35.67 ± 1.56 

34.02 ± 1.61 

35.18 

35.52 

35.22 

63.22 ± 1.60 

61.63 ± 1.65 

62.56 

63.43 

62.64 


B 

34.53 ± 1.13 

34.48 ± 1.22 

35.18 

35.39 


61.88± 1.17 

61.64 ± 1.25 

62.55 

63.11 



c 

35.05 ± 0.74 

35.05 ± 0.77 

35.15 

35.24 

±0.15 

62.39 ±0.76 

62.37 ±0.79 

62.59 

62.78 

±0.11 


D 

35.28 ± 0.54 

35.47 ± 0.53 

35.15 

35.19 


62.47 ±0.54 

62.97 ±0.55 

62.59 

62.68 


d 

A 

33.70 ± 2.02 

32.43 ± 2.14 

32.58 

32.76 

32.63 

61.44 ± 2.06 

63.26 ± 2.31 

62.84 

63.99 

62.97 

B 

31.45 ± 1.43 

32.26 ± 1.64 

32.58 

32.77 


62.58 ± 1.61 

62.41 ± 1.72 

62.90 

63.62 


-t 

c 

32.63 ± 0.96 

32.94 ± 0.99 

32.54 

32.68 

±0.19 

63.53 ± 1.02 

62.94 ± 1.06 

62.88 

63.22 

±0.20 


D 

32.31 ± 0.69 

33.00 ± 0.72 

32.52 

32.65 


62.55 ± 0.83 

62.43 ±0.86 

62.89 

63.09 




HMC 

SMC 

HPDE 

APDE 

A 

30 s 

30 s 

32 s 

30 s 

B 

122 s 

119 s 

131 s 

114 s 

C 

477 s 

476 s 

410 s 

491 s 

D 

1915 s 

1907 s 

1755 s 

1933 s 


2.5% 


2 . 0 % 


1.5% 

o 

W 

1 . 0 % 


0.5% 


0 . 0 % 


\ Test 2-Static: Relative Errors (error/time) 

\ 

A.... 



0 500 1000 1500 2000 Tinie (s) 
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—■*— HPDE 
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Table 7: Test 2-Static.In the first Table, the fair fee for the Heston model, with no ratchet or annual ratchet. 
In the second Table the run times for the no-ratchet case {p = —0.5). Finally, the plot of relative error (w.r.t. 
MC value) for the three methods in the case p = —0.5 with no ratchet. The parameters used for this test 
are available in Table [2] and in Table [H 
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Table 8: The model parameters about Test 2B-Static. 



No ratchet 

Annual ratchet 


HMC 

SMC 

HPDE 

APDE 

BM 

HMC 

SMC 

HPDE 

APDE 

BM 

A 

138.54± 2.70 

141.86 ± 2.70 

130.83 

137.13 

131.11 

125.40 ±2.46 

128.33± 2.77 

117.08 

124.19 

117.56 

B 

132.57± 1.98 

137.16± 2.18 

130.80 

135.76 


119.33 ± 1.80 

123.68± 2.00 

117.18 

122.78 


C 

131.10± 1.28 

135.79± 1.32 

130.80 

133.85 

±0.80 

117.74± 1.17 

124.49± 1.20 

117.23 

120.09 

±0.71 

D 

130.22 ±0.92 

132.17±0.90 

130.82 

133.02 


116.98 ±0.84 

118.72 ±0.82 

117.19 

119.62 



Table 9: Test 2B-Static. Fair fee for the Heston model, with no ratchet or annual ratchet. The parameters 
used for this test are available in Table |3] and in Table [8l 



p 


no ratchet (cig = 

i5o; 


ratchet (Ug = 250^ 

HMC 

SMC 


HPDE 

APDE 

MC 

HMC 

SMC 


HPDE 

APDE 

MC 



A 

6055 ± 12 

6060 ± 

8 

6058 

6034 


7123 ± 12 

7119 ± 

8 

7118 

7128 


-2 








6059 






7121 



B 

6043 ± 9 

6054 ± 

4 

6058 

6050 


7108 ± 9 

7115 ± 

4 

7119 

7126 




c 

6059 ± 6 

6057 ± 

2 

6058 

6052 


7117 ± 6 

7119 ± 

2 

7120 

7120 










±1 






±1 



D 

6059 ± 4 

6057 ± 

1 

6058 

6055 


7120 ± 4 

7119 ± 

1 

7120 

7120 




A 

6057 ± 13 

6057 ± 

9 

6060 

6026 


7390 ± 13 

7392 ± 

9 

7380 

7394 










6059 






7393 

"o 


B 

6052 ± 10 

6057 ± 

5 

6059 

6044 


7382 ± 10 

7389 ± 

5 

7387 

7391 



o 














o 


c 

6058 ± 6 

6057 ± 

2 

6058 

6050 


7389 ± 6 

7392 ± 

2 

7390 

7389 










±1 






±1 

o 


D 

6059 ± 5 

6058 ± 

1 

6058 

6055 


7391 ± 5 

7392 ± 

1 

7391 

7390 




A 

6095 ± 13 

6093 ± 

9 

6100 

6002 


7647 ± 14 

7651 ± 

9 

7636 

7649 


CQ 








6097 






7650 



B 

6101 ± 10 

6097 ± 

5 

6098 

6041 


7646 ± 10 

7648 ± 

5 

7643 

7644 



± 

c 

6098 ± 7 

6096 ± 

2 

6097 

6063 


7647 ± 7 

7651 ± 

2 

7647 

7644 










±1 






±2 



D 

6099 ± 5 

6097 ± 

1 

6097 

6080 


7650 ± 5 

7650 ± 

1 

7649 

7646 







no ratchet (ug = 

5o; 


ratchet (ag = 100^ 


P 


HMC 

SMC 

HPDE 

APDE 

MC 

HMC 

SMC 

HPDE 

APDE 

MC 


lO 

A 

7870 ± 20 

7856 ± 23 

7875 

7867 

7875 

8509 ± 14 

8499 ± 15 

8502 

8512 

8509 


B 

7873 ± 15 

7868 ± 16 

7875 

7873 


8506 ± 10 

8503 ± 11 

8506 

8516 




c 

7874 ± 9 

7877 ± 10 

7875 

7874 

±1 

.8505 ± 7 

8511 ± 7 

8507 

8512 

±1 



D 

7888 ± 7 

7880 ± 7 

7875 

7872 


8513 ± 5 

8513 ± 5 

8508 

8506 


O 

O! 

0) 

K 


A 

7803 ± 23 

7181 ± 25 

7797 

7786 

7897 

8405 ± 16 

8390 ± 17 

8395 

8400 

8398 


B 

7792 ± 16 

7790 ± 18 

7797 

7794 


8398 ± 12 

8391 ± 13 

8397 

8405 




c 

7796 ±11 

7801 ± 11 

7797 

7797 

±2 

8399 ± 8 

8398 ± 8 

8397 

8401 

±2 



D 

7803 ± 8 

7803 ± 8 

7797 

7795 


8402 ± 6 

8403 ± 6 

8398 

8395 



lO 

A 

7730 ± 31 

7719 ± 31 

7718 

7699 

7718 

8268 ± 22 

8292 ± 22 

8281 

8283 

8282 


B 

7703 ± 20 

7717 ± 22 

7717 

7712 


8292 ± 16 

8281 ± 16 

8282 

8290 



± 

c 

7718 ± 13 

7726 ± 14 

7717 

7717 

±3 

8283 ± 10 

8484 ± 10 

8282 

8286 

±2 



D 

7714 ± 9 

7723 ± 11 

7717 

7715 


8278 ± 7 

8287 ± 7 

8282 

8279 



Table 10: Test 3-Static. Delta calculation for the Black-Scholes Hull-White model and the Heston model, 
with no ratchet or annual ratchet (these value must be multiplied by 10“"^). The parameters used for this 
test are available in Table in Table and in Table 
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No ratchet 

Ratchet 



HMC 

SMC 

HPDE 

APDE 

BM 

HMC 

SMC 

HPDE 

APDE 

BM 



41.89 

41.75 

41.75 

41.85 

41.84 

75.96 

75.91 

75.92 

76.20 

76.00 

K 

-0.5 

±0.37 

±0.08 



±0.12 

±0.42 

±0.09 



±0.13 

0 

76.73 

76.21 

76.18 

76.22 

76.28 

143.12 

143.03 

142.95 

143.23 

143.06 

CO 

±0.48 

±0.11 



±0.15 

±0.56 

±0.12 



±0.18 



104.65 

104.48 

104.42 

103.89 

104.41 

204.49 

204.45 

204.24 

204.49 

204.37 


+0.5 

±0.61 

±0.13 



±0.19 

±0.73 

±0.16 



±0.23 




No ratchet 

Ratchet 



HMC 

SMC 

HPDE 

APDE 

BM 

HMC 

SMC 

HPDE 

APDE 

BM 


-0.5 

34.48 

34.19 

33.88 

33.92 

33.88 

56.03 

56.83 

55.54 

55.61 

55.56 



±0.43 

±0.42 



±0.11 

±0.42 

±0.41 



±0.11 

0 

0 

32.03 

32.23 

31.90 

31.94 

31.88 

55.89 

55.18 

55.80 

55.87 

55.85 

CD 

K 


±0.53 

±0.53 



±0.14 

±0.55 

±0.54 



±0.15 


±0.5 

28.97 

29.67 

29.25 

29.32 

29.26 

55.11 

55.93 

55.40 

55.53 

55.39 



±0.68 

±0.71 



±0.19 

±0.72 

±0.72 



±0.20 


Table 11: Test 4-Static. Impact of -1-10% mortality shocks of fair fee. The parameters used for this test are 
available in Table in Table and in Table 
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5.2 Dynamic case 


In the Dynamic case, the policy holder is supposed to choose the worst strategy from an hedger point of 
view, changing the value of 7^. The PH can withdraw more (1 < 7t < 2) or less (0 < 7t < 1) than the 
standard rate (see 2.4.3 for more details). 

In this pricing framework, we refer to the prices in Forsyth and Vetzal [H]: in their article, the authors 
price a GLWB contract in a static framework, under the Black Scholes model with r = 0.04 and a = 0.15. 
The contract parameters are reported in the Table 12 (Table 6.7 in j9j). They treated two cases: no ratchet, 
and ratchet every 3 years; both of them corresponds to case 4 in Section |4.1.3| and |4.3.4| In the first case 
they got ttg = 63.1 bp and in the second case ag = 70.7 bp. In Test 1 and 2 we introduce respectively 
stochastic interest rate and stochastic volatility to analyze the impact of these model developments on the 
fair guarantee fee. The parameters for interest rate and volatility models are the same as the static case. 

Here a brief summary of the numerical results for this Section. 


5.2.1 Test 1-Dynamic: the Black-Scholes Hull-White Model 

Test 1-Dynamic is the dynamic case of Test 1-Static. Model parameters are shown in Table Results are 
available in Table [13] 

In this Test, PDF methods proved to be more efficient than MC ones. In fact MC ones use Longstaff- 
Schwartz method to find the optimal withdrawal: this method needs a lot of scenarios to approximate 
through the least squares approach the value of the police for a given set of variable, and the regression is 
time demanding. Then, working at fixed time, we could perform fewer scenarios than static case (around 
10%), while PDF methods used almost the same parameters as in static case. Moreover the regression 
problem proved to be be hard: sometimes, excluding the value 7 = 0 among the possible values that the 
PH can chose (therefore excluding no withdrawal case), we got higher values for ag. This means that the 
regression isn’t very accurate, and sometimes we fail to find the optimal withdrawal: that’s why, using MC 
methods we usually find smaller value for ag than the right value. In particular, we excluded the value 
7 = 0 while using configurations A and B. We would remark that also the benchmark is affected by these 
computation problems and in case no-ratchet with p = —0.5 we got a small value for benchmark then PDF 
method (around 261 vs 266). Another thing to remark is that MC methods behaved better while ratchets 
were considered: maybe in this case in it easier to find the best strategy. The two MC methods proved to be 
equivalent: the differences in scenario generation’s run-time are negligible because most of the time is spent 
in finding the best withdrawal. Both APDF and HPDF method gave good and stable results, but HPDF 
performed better in case A. 


5.2.2 Test 2-Dynamic: the Heston Model 

Test 2-Dynamic is the dynamic case of Test 2-Static. Model parameters are shown in Table [^ Results are 
available in Table dn 

In this Test, things are similar to Test 1-Dynamic, but the optimization problem seemed to be easier than 
in Test 1-Dynamic: MC methods converged better, especially when using high level configurations. PDF 
methods behaved good as usual, and HPDF method proved to be a bit better then APDF method. The 


Initial age of PH 

65 

Gr. premium 

100 

DB payment 

cont.ly 

Strategy 

Dynamic 

G 

0.05 

Initial fees 

0 

Ratchet 

Off/On 

Bonus 

5% 

Withdrawal rate 

First withdrawal 

1 per Y 

P* anniv. 

Fees taken 

0 

cont.ly 

Ratchet rate 

every 3 Ys 

K{t) 

see tab below 


K{t) 

0 < t < 1 

1 < t < 2 

2 < t < 3 

3 < t < 4 

4 < t < 5 

t > 5 

5% 

4% 

3% 

2% 

1% 

0% 


Table 12: The contract parameters used in the dynamic case. 
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two MC methods proved to be equivalent. We note that, in Heston model case, Dynamic strategy increase 
the value of ag less than in BS HW case: probably, playing on interest rate, let the PH to gain more than 
playing on volatility. 

5.2.3 Test 3-Dynamic: Hedging 

Test 3-Dynamic is the dynamic case of Test 3-Static. Results are available in Table [T5| 

In this Test, we got good results with all methods, but MC methods proved to be inaccurate while using 
configurations A and B. The range of possible values for Delta increased with regard to Test 3-Static. The 
two MC methods proved to be equivalent. 

5.2.4 Test 4-Dynamic: Risk Management 

Test 4-Dynamic is the dynamic case of Test 4-Static. Results are available in Table [T6| 

In this Test, we got similar results with regard to Test 4-Static: the fees reduced a little (around 20 bp 
in the BS HW model case and around 6 bp in the Heston model case). 

In Figure |5.1[ we present, as an example, the optimal strategy at time t = 1 in two different cases. We 
can see the best strategy at time t = 1. We can see how it is worth to lapse when the account value reaches 
high values, and especially when interest rate is high or volatility low. It’s more difficult to understand when 
do no withdrawal: there must be a convenient mix of all the variables. 
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no ratchet 

ratchet 

p 


HMC 

SMC 

HPDE 

APDE 

BM 

HMC 

SMC 

HPDE 

APDE 

BM 

ic 

A 

79.71 ± 9.68 

66.5 ± 10.1 

85.50 

81.92 

84.35 

98.5 ± 16.1 

90.24 ± 10.5 

96.57 

92.80 

96.37 

B 

77.21 ± 5.10 

72.78 ± 4.52 

85.34 

84.78 


92.31 ± 5.12 

92.69 ± 5.01 

96.34 

95.75 



c 

79.86 ± 2.73 

80.29 ± 2.60 

85.27 

84.86 

±0.66 

94.70 ± 2.77 

93.79 ± 2.67 

96.25 

95.78 

±0.64 


D 

81.66 ± 1.58 

81.58 ± 1.46 

85.23 

84.54 


93.66 ± 1.62 

94.78 ± 1.62 

96.19 

95.40 



A 

162.6 ± 18.5 

148.4 ± 13.2 

172.55 

167.86 

169.05 

182.3 ± 13.2 

179.5 ± 14.3 

186.44 

181.96 

186.53 


B 

155.43 ±7.70 

156.52 ± 6.95 

172.60 

171.48 


182.42 ±6.16 

181.35±6.53 

186.48 

185.33 



c 

161.53 ±4.39 

164.88±4.41 

172.57 

171.61 

±0.90 

184.21 ±3.70 

183.84 ±3.58 

186.54 

185.45 

±0.86 


D 

164.51 ±2.53 

163.15±2.23 

172.58 

171.15 


183.87±2.21 

183.27±2.15 

186.55 

185.00 


ic 

A 

256.6 ± 15.6 

238.6 ± 23.3 

265.33 

261.22 

261.29 

273.9 ± 24.8 

262.7 ± 24.9 

272.18 

268.10 

274.02 

B 

246.0 ± 10.4 

248.5 ± 10.9 

266.66 

265.71 


268.37 ±8.39 

269.82±9.11 

273.24 

272.33 


-1- 

c 

D 

253.70 ±5.96 

253.33± 5.38 

266.93 

265.94 

±1.23 

271.90±5.51 

271.60 ±4.49 

273.67 

272.46 

±1.20 


259.00 ±3.38 

254.06± 3.11 

267.29 

265.38 


272.24± 2.95 

270.35 ±2.86 

273.99 

271.94 




HMC 

SMC 

HPDE 

APDE 

A 

30 s 

31 s 

32 s 

30 s 

B 

119 s 

122 s 

127 s 

120 s 

C 

482 s 

487 s 

463 s 

466 s 

D 

1911 s 

1942 s 

1732 s 

1815 s 


12 . 0 % 

10 . 0 % 

8 . 0 % 


Test 1-Dynamic: Relative Errors (error/time) _ 
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Table 13: Test 1-Dynamic. In the first Table, the fair fee for the Black-Scholes Hull-White model, with no 
ratchet or annual ratchet. In the second Table the run times for the no-ratchet case {p = —0.5). Finally, 
the plot of relative error (w.r.t. MC value) for the three methods in the case p = —0.5 with no ratchet. The 
parameters used for this test are available in Table 1^ and in Table 
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no ratchet 

ratchet 

p 


HMC 

SMC 

HPDE 

APDE 

BM 

HMC 

SMC 

HPDE 

APDE 

BM 

ic 

A 

56.96 ± 6.01 

62.56 ± 8.78 

64.57 

64.89 

65.38 

66.72 ± 8.37 

72.06 ± 6.17 

71.23 

71.58 

72.03 

B 

61.68 ± 3.05 

57.20 ± 3.96 

64.72 

64.64 


72.53 ± 3.56 

68.45 ± 3.36 

71.37 

71.30 



c 

64.68 ± 2.02 

64.03 ± 2.05 

64.76 

64.42 

±0.45 

71.65 ± 1.92 

71.95 ± 1.99 

71.43 

71.04 

±0.45 


D 

63.85 ± 1.28 

64.67 ± 1.28 

64.81 

64.35 


70.39 ± 1.22 

71.33 ± 1.25 

71.50 

70.96 



A 

58.83 ± 11.0 

58.68 ± 19.68 

61.92 

61.92 

62.32 

65.95 ± 9.04 

79.17 ± 10.36 

68.97 

68.98 

69.54 


B 

58.95 ± 4.01 

54.34 ±4.08 

61.91 

61.67 


69.14 ±4.53 

65.30 ±4.25 

68.95 

68.69 



c 

58.26 ± 2.25 

57.76 ± 2.34 

61.88 

61.43 

±0.56 

68.89 ± 2.43 

68.50 ± 2.44 

68.94 

68.41 

±0.58 


D 

59.16 ± 1.43 

59.70 ± 1.47 

61.87 

61.35 


66.76 ± 1.53 

68.67 ± 1.57 

68.93 

68.33 


ic 

A 

52.66 ± 15.29 

61.70 ± 13.04 

57.25 

57.50 

56.59 

63.67 ± 9.76 

83.75 ± 13.01 

64.60 

64.79 

65.42 

B 

57.26 ± 5.06 

51.95 ± 6.09 

57.33 

57.26 


68.21 ± 6.30 

60.31 ± 5.94 

64.66 

64.53 


-1- 

c 

52.23 ± 3.01 

51.47 ±3.50 

57.36 

57.01 

±0.73 

63.78 ± 3.13 

64.06 ± 3.25 

64.68 

64.25 

±0.74 


D 

52.60 ± 1.85 

52.24 ± 1.83 

57.39 

56.94 


62.98 ± 1.87 

61.86 ± 1.97 

64.71 

64.16 



. lest 2-Dynamic: Helative Errors (error/time) 
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Table 14: Test 2-Dynamic. In the first Table, the fair fee for the Heston model, with no ratchet or annual 
ratchet. In the second Table the run times for the no-ratchet case {p = —0.5). Finally, the plot of relative 
error (w.r.t. MC value) for the three methods in the case p = —0.5 with no ratchet. The parameters used 
for this test are available in Table [12] and in Table [6| 
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no ratchet (ctg = 300) 


ratchet (a 

3 = 350) 




P 


HMC 

SMC 

HPDE 

APDE 

BM 

HMC 

SMC 

HPDE 

APDE 

BM 



A 

8513 ± 573 

7951 ± 557 

8078 

8097 


8959 ± 542 

8279 ± 551 

8148 

8200 









8078 





8157 



B 

8220 ± 347 

8304 ± 305 

8081 

8091 


7922 ± 324 

8104 ± 330 

8151 

8180 




c 

8091 ± 172 

8146 ± 161 

8082 

8073 


8155 ± 173 

8120 ± 162 

8152 

8163 









±42 





±38 

=! 


D 

8089 ± 108 

8105 ± 104 

8082 

8093 


8164 ± 101 

8174 ± 99 

8152 

8164 




A 

7898 ± 550 

7697 ± 489 

7516 

7531 


7733 ± 509 

7379 ± 667 

7538 

7539 









7485 





7517 

"o 


B 

7685 ± 257 

7389 ± 269 

7517 

7529 


7640 ± 259 

7417 ± 249 

7527 

7539 



o 












o 


c 

7488 ± 139 

7443 ± 150 

7517 

7514 


7433 ± 145 

7604 ± 137 

7528 

7530 









±31 





±29 

o 


D 

7428 ± 90 

7489 ± 83 

7517 

7518 


7523 ± 87 

7525 ± 84 

7528 

7533 




A 

7444 ± 470 

7612 ± 491 

7333 

7342 


7569 ± 421 

7292 ± 500 

7304 

7314 


CQ 







7324 





7309 



B 

7257 ± 242 

7368 ± 209 

7337 

7350 


7386 ± 207 

7413 ± 192 

7308 

7322 



+ 

c 

7306 ± 116 

7201 ± 124 

7339 

7336 


7469 ± 116 

7293 ± 112 

7309 

7307 









±28 





±26 



D 

7270 ± 78 

7302 ± 78 

7340 

7339 


7291 ± 70 

7298 ± 68 

7310 

7310 





no ratchet (ag = 75) 


ratchet (a 

9 = 100) 



P 


HMC 

SMC 

HPDE 

APDE 

BM 

HMC 

SMC 

HPDE 

APDE 

BM 

lO 

A 

8181 ± 472 

7794 ± 524 

8436 

8429 

8432 

8349 ± 293 

8374 ± 324 

8477 

8474 

8481 

B 

8440 ± 250 

8383 ± 223 

8436 

8436 


8304 ± 174 

8527 ± 154 

8479 

8480 



c 

8405 ± 87 

8426 ± 87 

8437 

8437 

±19 

8535 ± 80 

8499 ± 76 

8479 

8480 

±14 


D 

8437 ± 50 

8472 ± 53 

8437 

8438 


8516 ± 44 

8501 ± 43 

8480 

8480 



A 

8756 ± 626 

8304 ± 586 

8329 

8319 

8297 

8751 ± 758 

8440 ± 420 

8341 

8332 

8351 


B 

8080 ± 283 

8345 ± 394 

8330 

8327 


8466 ± 251 

8184 ± 217 

8341 

8339 



c 

8313 ± 148 

8137 ± 185 

8330 

8329 

±29 

8303 ± 106 

8398 ± 101 

8341 

8340 

±18 


D 

8330 ± 75 

8238 ± 79 

8330 

8331 


8343 ± 52 

8283 ± 66 

8341 

8441 


lO 

A 

7308 ± 1145 

7453 ± 914 

8217 

8205 

8242 

8244 ± 736 

8192 ± 522 

8191 

8180 

8206 

B 

8238 ±623 

7919 ±416 

8218 

8215 


8150 ± 276 

8213 ± 229 

8191 

8189 


± 

c 

8143 ± 145 

7874 ± 241 

8218 

8217 

±46 

8216 ± 133 

8123 ± 178 

8192 

8190 

±22 


D 

8144 ± 119 

8131 ± 108 

8218 

8219 


8123 ± 66 

8195 ± 77 

8192 

8191 



Table 15: Test 3-Dynamic. Delta calculation for the Black-Scholes Hull-White model and the Heston model, 
with no ratchet or ratchet (once every 3 years); these value must be multiplied by 10“"^. The parameters 
used for this test are available in Table in Table and in Table 
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No ratchet 

Ratchet 



HMC 

SMC 

HPDE 

APDE 

BM 

HMC 

SMC 

HPDE 

APDE 

BM 



76.64 

75.74 

76.85 

76.19 

76.16 

85.10 

85.47 

86.59 

85.84 

86.97 

& 

s 

-0.5 

±1.49 

±1.37 



±0.60 

±1.53 

±1.51 



±0.61 

0 

149.80 

149.34 

155.47 

154.11 

155.97 

168.14 

167.94 

169.78 

168.26 

170.32 


±2.38 

±2.14 



±0.96 

±2.14 

±2.00 



±0.96 



236.49 

236.33 

242.08 

240.22 

242.1 

248.40 

246.49 

250.08 

248.12 

250.24 


+0.5 

±2.84 

±2.79 



±1.2 

±2.68 

±2.48 



±0.53 




No ratchet 

Ratchet 



HMC 

SMC 

HPDE 

APDE 

BM 

HMC 

SMC 

HPDE 

APDE 

BM 



57.74 

58.77 

58.78 

58.32 

58.56 

63.73 

64.50 

64.74 

64.22 

64.59 


-0.5 

±1.26 

±1.27 



±0.46 

±1.23 

±1.27 



±0.46 

0 

0 

53.08 

53.77 

55.57 

55.07 

55.76 

59.86 

61.17 

61.79 

61.22 

62.11 

0) 

±1.45 

±1.46 



±0.58 

±1.51 

±1.59 



±0.59 



46.33 

46.19 

50.93 

50.50 

50.66 

55.67 

54.79 

57.29 

56.76 

59.01 


+0.5 

±1.87 

±1.77 



±0.70 

±1.85 

±1.94 



±0.76 


Table 16: Test 4-Dynamic. Impact of -1-10% mortality shocks of fair fee. The parameters used for this test 
are available in Table in Table and in Table 
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Figure 5.1: Optimal strategy at the first event time {t = 1) for the BS HW model and the Heston model, 
assuming = 100. Model parameters are available in Tables and Product parameters are available 


in Table 12 and ag = 135 bp for both cases. 
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6 Conclusions 


In this article we have developed four methods to price GLWB contracts under different conditions. Regarding 
the stochastic model, both stochastic interest rates and stochastic volatility effects have been considered. 
Regarding the policy holder behavior, both static and dynamic strategy have been considered. 

Since GLWB variable annuities are such a long maturity products, the effects of stochastic interest rates 
and stochastic volatility cannot be overlook. In particular, the impact of stochastic rate seems to be more 
relevant. Also Forsyth and Vetzal in [5] used regime switching model having both stochastic interest rate 
and volatility, but our approach, based on SDE, is more realistic, and suitable for hedging. 

All four methods gave compatible results both for pricing and delta calculation. The fair hedging fee (i.e. 
the cost of maintaining the replicating portfolio) is determined using a sequence of parameters’ refinements. 
The PDF methods proved to be not very expensive, while MG methods proved to be more expensive. The 
Hybrid PDF seemed to be the more performing than the others for its convergence speed and stability of 
results. Also ADI PDF behaved very well but the implementation was harder then Hybrid PDF one. In the 
BS HW model case. Standard MG thanks to its exact simulation outperformed the Hybrid Method while, in 
the Heston model case, the MG methods proved to be roughly equivalent, even if the Hybrid MG was easier 
to be implemented. 

As we said before, PDF methods proved to be much more efficient than MG methods, especially in 
Dynamic case where is much more simple to implement the optimal withdrawal choice. Similarity reduction 
reduces the problem dimension to a 2D problem and therefore PDF methods perform well. Anyway, we 
have to remark that MG methods offer a confidence interval for the result, they are useful in risk measures 
calculation (for example VAR or FS), and they are preferred by insurance companies because of their 
attachment to the idea of scenario. 

A future development that could be treated is to combine stochastic interest rate and stochastic volatility: 
the combined model could be an element of greater realism. 

We conclude by pointing out that our methods are quite flexible in that they can accommodate a wide 
variety of policy holder withdrawal strategies such as ones derived from utility-based models. 


32 


References 


[1] A. Alfonsi (2010). High order discretization schemes for the CIR process: application to Affine Term Structure 
and Heston models. Mathematics of Computation, Vol. 79, No. 269, pp. 209-237. 

[2] E. Appolloni, L. Caramellino A. Zanette 2014. A robust tree method for pricing American options with 
the Cox-Ingersoll-Ross interest rate model. IMA J Management Math first published online January 15, 2014 
doi: 10.1093/imaman / dpt030. 

[3] A. R. Bacinello, P. Millossovich, A. Olivieri, E. Pitacco (2011). Variable annuities: A unifying valuation 
approach. Insurance: Mathematics and Economics 49, pp. 285-297. 

[4] M. Briani, L. Caramellino, A. Zanette (2014). A hybrid tree-finite difference approach for the Heston 
model. Preprint, arXiv: 1307.7178 i^2. 

[5] M. Briani, L. Caramellino, A. Zanette (2015). Numerical approximations for Heston-Hull-White type 
models. Preprint, arXiv:1503.03705 . 

[6] A. Belanger, P. Forsyth, G. Labahn (2009). Valuing the guaranteed minimum death benet clause with 
partial withdrawals. Applied Mathematical Finance 16, pp. 451-496. 

[7] D. Brigo, F. Mercurio (2006). Interest rate models-Theory and practice. Springer, Berlin. 

[8] Z. Chen, K. Vetzal, P. Forsyth (2008). The effect of modelling parameters on the value of GMWB guarantees. 
Insurance: Mathematics and Economics 43, pp. 165-173. 

[9] P. Forsyth, K. Vetzal (2014). An optimal stochastic control framework for determining the cost of hedging of 
variable annuities. Journal of Economic Dynamics and Control 44 (2014), pp. 29-53. 

[10] P. Gaillardetz, J. Lakhmiri (2011). A new premium principle for equity indexed annuities. Journal of Risk 
and Insurance 78, 245-265. 

[11] S. Heston (1993): A closed-form solution for options with stochastic volatility with applications to bond and 
currency options. The Review of Financial Studies, Vol. 6, No. 2, pp. 327-343. 

[12] D. Holz, a. Kling, j. Russ (2007). GMWB for life: an analysis of lifelong withdrawal guarantees. Working 
paper. 

[13] T. Haentjens, K. J. In ’t Hout (2012): Alternating direction implicit finite difference schemes for the 
Heston-Hull-WHite partial differential equation. The Journal of Gomputation Finance (83-110), Vol. 16, No. 1, 
Fall 2012. 

[14] J. Hull, A. White (1994). Numerical procedures for implementing term structure models I: single factor 
models. The Journal of Derivatives Fall, 716. 

[15] A. Kling, F. Ruez, J. Russ (2014). The impact of stochastic volatility on pricing, hedging, and hedge efficiency 
of variable annuity guarantees. European Actuarial Journal, Vol. 4, No. 2, pp. 281-314. 

[16] F. A. Longstaff, E. S. Schwartz (2001). Valuing american options by simulation: a simple least-squares 
approach. The Review of Financial Studies Spring 2001, Vol. 14, No. 1, pp. 113-147. 

[17] M. A. Milevsky, T. S. Salisbury (2006). Financial valuation of guaranteed minimum withdrawal benets. 
Insurance: Mathematics and Economics 38, pp. 21-38. 

[18] D. B. Nelson, K. Ramaswamy (1990). Simple binomial processes as diffusion approximations in financial 
models. The Review of Financial Studies 1990, Vol. 3, No. 3, pp. 393-430. 

[19] P. Shah, D. Bertsimas (2008). An analysis of the guaranteed withdrawal benets for life option. Working paper, 
Sloan School of Management, MIT. 

[20] V. Ostrovski (2013). Efficient and exact simulation of the Hull-White model. Available at SSRN: 
http://ssrn.com/abstract=2304848 or http://dx.doi.org/10.2139/ssrn.2304848. 


33 



